import scvelo as scv
scv.logging.print_version()
scv.settings.verbosity = 3 # show errors(0), warnings(1), info(2), hints(3)
scv.settings.presenter_view = True # set max width size for presenter view
scv.set_figure_params('scvelo') # for beautified visualization
- The analysis is based on the in-built pancreas data.
- To run velocity analysis on your own data, read your file (loom, h5ad, csv …) to an AnnData object with
adata = scv.read('path/file.loom', cache=True). - If you want to merge your loom file into an already existing AnnData object, use scv.utils.merge(adata, adata_loom).
adata = scv.datasets.pancreas()
adata
Names of observations and variables can be accessed via adata.obs_names and adata.var_names, respectively. AnnData objects can be sliced like dataframes, for example, adata_subset = adata[:, list_of_gene_names].
scv.pl.proportions(adata)
- ここでは、スプライスされた/されていない数の割合が表示されます。使用したプロトコル（Drop-Seq、Smart-Seq）にもよりますが、通常、intron配列を含むスプライスされていない分子が10％～25％程度存在します。
- また、スプライシング効率の一貫性を確認するために、クラスターレベルでのばらつきを調べることをお勧めします。ここでは、予想通りのばらつきが見られ、周期性管細胞 (cycling ductal cells)ではスプライスされていない分子の割合がわずかに低く、その後、多くの遺伝子が転写され始めるNgn3-high細胞やPre-endocrine細胞では、細胞運命コミットメントでの割合が高くなっていることがわかりました。
Preprocess the Data¶
- 遺伝子選択(最小発現量とhigh variability (dispersion))、
- 細胞をTotal countで正規化し、対数化
フィルタリングと正規化は、スプライスされた/されていないカウントにも同じように適用されます。
これらはすべて単一の関数 scv.pp.filter_and_normalize にまとめられており、以下のコマンドが含まれます。
scv.pp.filter_and_normalize(adata, min_shared_counts=20, n_top_genes=2000)
#scv.pp.filter_genes(adata, min_shared_counts=20)
#scv.pp.normalize_per_cell(adata)
#scv.pp.filter_genes_dispersion(adata, n_top_genes=2000)
#scv.pp.log1p(adata)
- さらに、PCA空間の最近傍点間で計算された1次と2次のモーメント（平均と非中心分散）が必要で、
scv.pp.momentsにまとめられています (内部でscv.pp.pcaとscv.pp.neighborsを利用)。 - 決定論的速度推定には1次が必要ですが，確率的推定には2次モーメントが必要です．
scv.pp.moments(adata, n_pcs=30, n_neighbors=30)
さらなる前処理 (optional)¶
- さらなる前処理（バッチ効果補正など）を使用して、不要な変動源を除去することができます。
- 追加の前処理ステップはXにのみ影響し、スプライス/アンスプライス数には適用されないことに注意してください。
Estimate RNA velocity¶
- Velocityは遺伝子発現空間におけるベクトルであり、個々の細胞の移動の「方向」と「速度」を表します。
- Velocity = speed + direction
- Velocitiesは、確率的（デフォルト）または決定論的 (by setting mode='deterministic')に、スプライシング動態の転写ダイナミクスをモデル化することで得られます。
- 各遺伝子について、pre-mature (unspliced), mature (spliced)の mRNAカウントの定常状態比 (steady-state-ratio) をフィットさせ、これが一定の転写状態を表す。
- そして、Velocityはこのsteady-state-ratioの残差として求めます。
- 正のVelocityは、ある遺伝子が発現上昇していることを示しており、定常状態で予想されるよりも多くunspliced mRNAが存在する細胞で発生します。
逆に、負の速度は遺伝子がダウンレギュレーションされていることを示している。
完全な動的モデルの解は、
mode='dynamical'を設定することで得られます。このモードは事前にscv.tl.recover_dynamics(adata)を実行しておく必要があります。このモデルについては次のチュートリアルで詳しく説明します。
scv.tl.velocity(adata)
計算された速度は、カウント行列と同様に adata.layer に保存されます。
- 遺伝子間のvelocityの組み合わせは、個々の細胞の将来の状態を推定するために使用することができます。
- 速度を低次元の埋め込みに投影するために、細胞間遷移確率が推定されます。つまり、各速度ベクトルについて、その方向に沿ったセル遷移の可能性が高いものを見つけます。
- 遷移確率は、潜在的な細胞間遷移と速度ベクトルとの間の余弦相関を用いて計算され、velocity graphと呼ばれる行列に格納される。
- 結果として得られるvelocity graphは$n_{obs} \times n_{obs}$ の次元を持ち、速度ベクトルによって説明できる細胞の状態変化を要約しています（高速化のために、
approx=Trueを設定することで、縮小されたPCA空間上で計算することもできます）。
scv.tl.velocity_graph(adata)
- 様々なアプリケーションのために、速度グラフを遷移行列に変換することができます。これはガウスカーネルを適用してコサイン類似度を実際の遷移確率に変換することで行います。
マルコフ遷移行列には
scv.utils.get_transition_matrixでアクセスできます。前述したように、velocityを低次元空間に埋め込むために、
scv.tl.velocity_embeddingを使用して得られた遷移確率に対する平均遷移を計算します。- さらに、マルコフ連鎖に沿って細胞をその起源と潜在的な運命までトレースすることができ、それにより、軌跡の根と終点を得ることができます。これは
scv.tl.terminal_statesで参照できます。
Project the velocities¶
最後に、速度は任意の低次元軸に投影され、以下のいずれかの方法で可視化されます。
- on cellular level with
scv.pl.velocity_embedding, - as gridlines with
scv.pl.velocity_embedding_grid, - or as streamlines with
scv.pl.velocity_embedding_stream.
- on cellular level with
デモデータには既に計算済みの UMAP 埋め込みと注釈付きクラスタが含まれています。
- 独自のデータに適用する場合、これらは
scv.tl.umapとscv.tl.louvainで取得できます。詳細については、scanpyのチュートリアルを参照してください。 - さらに、すべてのプロット関数はデフォルトで
basis='umap'とcolor='cluster'を使用するように設定されています。
scv.pl.velocity_embedding_stream(adata, basis='umap')
- ストリームラインとして表示されるvelocity ベクトルフィールドは、発生過程の詳細な洞察をもたらします。
- 管腔細胞(ductal cells) と内分泌前駆細胞 (endocrine progenitors) の周期的集団を正確に描き出します。
- さらに、細胞の系統コミットメント、細胞周期の出口、および内分泌細胞の分化の細胞の状態を表す。
cell level velocity¶
- 速度ベクトル場の最も詳細な分解能は、単一細胞レベルで得られるもので、各矢印は個々の細胞の移動方向と速度を示している。
- これにより、例えば、Ngn3細胞の初期の内分泌関与（黄色）や、近末端のα細胞（青）と一過性のβ細胞（緑）の間の明確な違いが明らかになりました。
scv.pl.velocity_embedding(adata, arrow_length=3, arrow_size=2, dpi=120)
Interprete the velocities¶
これはおそらく最も重要な部分です。つまり、生物学的な結論を投影されたvelocityに限定するのではなく、推定された方向が特定の遺伝子にどのように支持されているかを理解するために、相図 (phase portraits)を介して個々の遺伝子のダイナミクスを調べるべきです。
spliced vs. unspliced phase portrait をどのように解釈するかについては、こちらのgif を参照してください。
- 遺伝子の活性は転写制御によって調整されている。
- 特定の遺伝子に対して転写が誘導されると、（新たに転写された）precursor unspliced mRNAsが増加し、逆に転写が抑制されたり、転写が行われなかったりすると、unspliced mRNAsが減少します。
- スプライスされたmRNAは、未転写mRNAから産生され、時間差を持って同じ傾向をたどる。
- 時間は、 hidden/latent 変数である。
- そのため、実際に測定されているもの、すなわち、相図に表示されているspliced and unspliced mRNAsからダイナミクスを推測する必要があります。
さて、いくつかのマーカー遺伝子の相図を scv.pl.velocity(adata, gene_names) または scv.pl.scatter(adata, gene_names) で可視化してみましょう。
scv.pl.velocity(adata, ['Cpe', 'Gnao1', 'Ins2', 'Adk'], ncols=2)
- 黒線は推定された ‘steady-state’ の比率、すなわち、一定の転写状態における unspliced to spliced mRNA abundanceの比率に対応しています。
- 特定の遺伝子のRNA velocity は、残差、すなわち観測がその定常状態線 (黒線)からどれだけ逸脱しているかとして決定されます。
- 正の速度は、ある遺伝子が発現上昇していることを示しており、定常状態で予想されるよりもunspliced mRNAが豊富な細胞で発生します。
逆に、負の速度は遺伝子が発現抑制されていることを示している。
例えば、Cpeは発現上昇したNgn3（黄色）からPre-endocrine（オレンジ）からβ細胞（緑）への方向性を説明し、Adkは発現抑制されたDuctal（濃い緑）からNgn3（黄色）から残りの内分泌細胞への方向性を説明している。
scv.pl.scatter(adata, 'Cpe',
color=['clusters', 'velocity'],
add_outline='Ngn3 high EP, Pre-endocrine, Beta')
Identify important genes¶
私たちは、結果として得られるベクターフィールドと推定される系統を説明するのに役立つ遺伝子を特定する体系的な方法を必要としています。そのために、どの遺伝子がクラスター特異的な微分速度発現を持ち、残りの集団と比較して有意に高いか低いかを検定するモジュールを提供します。
scv.tl.rank_velocity_genesモジュールは differential velocity t-test を実行し、各クラスタの遺伝子ランキングを出力します。- 遺伝子候補選択のテストを制限するための閾値を設定することができます (例:
min_corr)。
scv.tl.rank_velocity_genes(adata, groupby='clusters', min_corr=.3)
df = scv.DataFrame(adata.uns['rank_velocity_genes']['names'])
df.head()
kwargs = dict(frameon=False, size=10, linewidth=1.5,
add_outline='Ngn3 high EP, Pre-endocrine, Beta')
scv.pl.scatter(adata, df['Ngn3 high EP'][:5], ylabel='Ngn3 high EP', **kwargs)
scv.pl.scatter(adata, df['Pre-endocrine'][:5], ylabel='Pre-endocrine', **kwargs)
The genes Ptprs, Pclo, Pam, Abcc8, Gnas, for instance, support the directionality from Ngn3 high EP (yellow) to Pre-endocrine (orange) to Beta (green).
Velocities in cycling progenitors¶
RNA速度によって検出される細胞周期は、細胞周期スコア（位相マーカー遺伝子の平均発現量を標準化したスコア）によって生物学的に肯定されます。
scv.tl.score_genes_cell_cycle(adata)
scv.pl.scatter(adata, color_gradients=['S_score', 'G2M_score'], smooth=True, perc=[5, 95])
周期性のある管状細胞 (Ductal cells)については、SとG2Mの位相マーカーを介してスクリーニングすることができます。上のモジュールはまた、位相マーカー遺伝子のランク付けやソートに使用できるスピアマン相関スコアを計算しており、それらのphase portraitsを表示するために使用することができます。
s_genes, g2m_genes = scv.utils.get_phase_marker_genes(adata)
s_genes = scv.get_df(adata[:, s_genes], 'spearmans_score', sort_values=True).index
g2m_genes = scv.get_df(adata[:, g2m_genes], 'spearmans_score', sort_values=True).index
kwargs = dict(frameon=False, ylabel='cell cycle genes')
scv.pl.scatter(adata, list(s_genes[:2]) + list(g2m_genes[:3]), **kwargs)
特にHellsとTop2aは、循環する前駆体のベクトル場を説明するのに適している。Top2aは、G2M段階で実際にピークを迎える直前に高い速度が割り当てられる。そこでは、負の速度は、その後のダウンレギュレーションの直後に完全に一致する。
scv.pl.velocity(adata, ['Hells', 'Top2a'], ncols=2, add_outline=True)
Speed and coherence¶
さらに2つの有用な統計量があります。
- 微分の速度や速度は速度ベクトルの長さで与えられます。
- ベクトル場のコヒーレンス（すなわち、速度ベクトルが隣接する速度とどのように相関するか）は、信頼性の尺度となります。
scv.tl.velocity_confidence(adata)
keys = 'velocity_length', 'velocity_confidence'
scv.pl.scatter(adata, c=keys, cmap='coolwarm', perc=[5, 95])
これらの結果から、細胞の分化のペースが遅い場合と早い場合、またその方向が不明な場合があることが明らかになった。
クラスターレベルでは、細胞周期の終了後（Ngn3低EP）には分化が実質的に加速され、ベータ細胞の産生中にはペースを維持し、アルファ細胞の産生中にはペースが低下することがわかった。
df = adata.obs.groupby('clusters')[keys].mean().T
df.style.background_gradient(cmap='coolwarm', axis=1)
Velocity graph and pseudotime¶
- Velocityグラフを可視化して、推定された全ての細胞間接続/遷移を描写することができます。
- 閾値を設定することで、高確率の遷移に限定することができる。
- 例えば、このグラフは、初期と後期の前内分泌細胞からのイプシロン細胞生産の2つのフェーズを示しています。
scv.pl.velocity_graph(adata, threshold=.1)
さらに、このグラフは、指定された細胞の子孫/先祖を描くために使用することができます。ここでは、前内分泌細胞は、その潜在的な運命までトレースされています。
x, y = scv.utils.get_cell_transitions(adata, basis='umap', starting_cell=70)
ax = scv.pl.velocity_graph(adata, c='lightgrey', edge_width=.05, show=False)
ax = scv.pl.scatter(adata, x=x, y=y, s=120, c='ascending', cmap='gnuplot', ax=ax)
- 最後に、速度グラフに基づいて、 velocity pseudotimeを計算することができる。
グラフからroot細胞からの分布を推定した後、root細胞からグラフに沿って歩いた後、その細胞に到達するまでにかかる平均ステップ数を測定することで計算されます。
diffusion pseudotimeとは逆に、root細胞を暗黙的に推定し、similarity-based diffusion kernelの代わりに有向速度グラフに基づいています。
scv.tl.velocity_pseudotime(adata)
scv.pl.scatter(adata, color='velocity_pseudotime', cmap='gnuplot')
PAGA velocity graph¶
- PAGAグラフ抽象化は、軌跡推論のための最高性能の手法としてベンチマーキングされています。
- PAGAは、2つのクラスタ間の接続性に対応した重み付きエッジを用いて、データのトポロジーをグラフ化したようなマップを提供する。
- ここでは、PAGAを速度推定の方向性によって拡張している。
scv.tl.paga(adata, groups='clusters')
df = scv.get_df(adata, 'paga/transitions_confidence', precision=2).T
df.style.background_gradient(cmap='Blues').format('{:.2g}')
これは、左/行から右/列へと読み取ることができ、例えば、DuctalからNgn3 low EPへの信頼性ある遷移を割り当てることができます。
この表は、UMAP上に埋め込まれた有向グラフでまとめることができます。
scv.pl.paga(adata, basis='umap',
size=50,
alpha=.1,
min_edge_width=2,
node_size_scale=1.5)
